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ABSTRACT 

We present a fitting formula for tlie non-linear evolution of the bispectrum in 
CDM models, obtained from measurements in high resolution numerical simulations. 
The formula interpolates between the perturbative and highly non-linear regimes, and 
generalizes previous results obtained for scale- free initial conditions. 

Key words: cosmology: large-scale structure of the universe; methods: numerical 



> 

O 

o 
o 

Oh' 

I ■ 



% 



1 INTRODUCTION 

The large-scale structure of the universe observed in galaxy 
surveys is thought to arise from the growth of primordial 
fluctuations due to gravitational instability. Making accu- 
rate predictions of gravitational clustering is essential if we 
are to compare theories with observations and put con- 
straints on cosmological parameters. 

Gravitational clustering is a non-linear process: when 
fluctuations are sufficiently small the equations of motion 
can be solved perturbatively, but at small enough scales one 
has resort to numerical simulations to accurately follow the 
dynamics. Non-linearity of the equations of motion implies 
that the mass distribution becomes non-Gaussian, even if 
the initial conditions are Gaussian (Peebles i980. Fry i984). 

The power spectrum, or two-point correlation function, 
has become the most common statistical measure of gravi- 
tational clustering. Its importance rests on the fact that for 
a Gaussian random field, the power spectrum completely 
characterizes its statistical properties. Hamilton et al. (I99I) 
suggested a useful way of thinking about the non- linear evo- 
lution of the power spectrum, summarized by a simple fitting 
formula obtained from numerical simulations. The accuracy 
of this procedure has been improved (e.g. Jain, Mo & White 
1995; Peacock & Dodds 1996) and extensively used in the 
literature to predict statistical properties related to the den- 
sity power spectrum. 

A complete description of the non- Gaussian distribu- 
tion that results from gravitational instability, however, re- 
quires knowledge of higher-order correlation functions. A 
number of useful statistics have emerged in the last decade 
to become standard probes of non-Gaussianity in large-scale 
structure. These include the higher-order moments of the 
smoothed density field (e.g. skewness and kurtosis; Peebles 
1980; Juszkiewicz, Colombi & Bouchet 1993; Bernardeau 



1994) and the three-point correlation function in Fourier 
space, the bispectrum (Fry 1984; Scoccimarro et al. 1998, 
hereafter SCFFHM) or in real space (Jing & Borner 1998; 
Frieman & Gaztahaga 1999). 

A fitting formula for the non-linear evolution of one- 
point higher-order moments is given by extended perturba- 
tion theory (EPT, Colombi et al. 1997). The EPT prescrip- 
tion uses the perturbative expressions for skewness, kurto- 
sis, and so on, but with an "effective" spectral index that 
becomes a function of the variance a^ of the smoothed den- 
sity field. This function is obtained numerically by fitting to 
measurements of the density field moments in N-body sim- 
ulations. The EPT prescription is able to fit reasonably well 
the non-linear evolution of higher-order moments and thus 
the one-point probability distribution function with this sin- 
gle function obtained from numerical simulations. 

Probing non-Gaussianity beyond one-point statistics 
gives precious additional information that can be used to 
learn about cosmology. The bispectrum is the lowest-order 
statistic that probes the shape of large-scale structures gen- 
erated by gravitational clustering. Measurements of the bis- 
pectrum in galaxy redshift catalogs have been proposed as 
a way to constrain galaxy biasing (Fry 1994; Hivon et al. 
1995; Matarrese, Verde & Heavens 1997) and thus break the 
degeneracies present in power spectrum measurements be- 
tween the matter density of the universe Q.m and galaxy bias. 
Such measurements have recently been done in IRAS red- 
shift surveys (Scoccimarro et al. 2000; Feldman et al. 2000), 
providing constraints on galaxy bias, Q,m and non-Gaussian 
primordial fiuctuations. In addition, upcoming large redshift 
surveys such as 2dF and SDSS will determine the bispectrum 
with unprecedented accuracy. 

In this work we study the description of the non-linear 
evolution of the bispectrum. A fitting formula for the bis- 
pectrum was proposed by Scoccimarro & Frieman (1999), 
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hereafter SF, that interpolates between perturbative and the 
highly non- linear regime, assuming scale- free initial condi- 
tions. Since current models for structure formation are not 
scale-free, these results could not be used to make realistic 
predictions to be compared with observations. In this work, 
we extend the SF results to provide a fitting formula for the 
non-linear evolution of the bispectrum in CDM models. 



2 THE BISPECTRUM 

We denote by S the Fourier transform of the density contrast 
S obtained from the density field p — p{l + S), where p is the 
mean density. The power spectrum P{k) and the bispectrum 
Bi23 = B{ki, fc2, fca) are defined as 



Table 1. Numerical Simulations Parameters 



(5(ki)5(k2)> 
(5(ki)5(k2)5(k3)> 



fc(ki-Fk2)P(fc) 

fc(ki + k2-f kg) -B123. (1) 



In this paper, we assume Gaussian initial conditions. In this 
case, the bispectrum generated by gravitational instability 
at large scales is given by tree-level (second-order) pertur- 
bation theory, and reads (Fry 1984) 



B123 = 2 F2(ki, ka) P1P2 + eye, 
where Pi = P{ki), and 



(2) 



5 lkiJC2/fcl^fc2\^2/kiJC2\2 

F2(k„k2) = - + -^^(- + -j+-(^^j .(3) 

Note that Eq.([l|) implies that the bispectrum is defined only 
for closed triangles formed by the wave vectors ki , k2 , ka . It 
is convenient to define the reduced bispectrum (3i23 as 

-B123 



Q123 



(4) 



P1P2 + P2P3 + P3P1 
which takes away most of the dependence on scale and cos- 
mology. 

SF gave a fitting formula for the non-linear evolution 
of the bispectrum from scale-free initial conditions. They 
replaced the kernel in Eq. (g) by the effective kernel 



F2 (ki,k2) = - a{n,ki)a{n,k2) 



+ 



7(^117) <^^^^^<^'^^)' 



(5) 



where the functions a{n,k), h{n,k), and c{n,k) are chosen 
to interpolate between these two regimes according to the 
one-loop FT and N-body results (SCFFHM). This yields 
(-2 < n < 0) 



a{n, k) — 
b{n, k) = 
c{n, k) — 



1 + [0.7 Q3(n)]^/^ {kRor+<' 
1 + (kRoy+^ 

1 + 0.2 {n + 2,) (fcJ^o)"+^ 
1 + (fc7?o)"+3-5 

l + ^S/iLS + ln + S)-*] (fci?o)"+' 



(6) 

(7) 
(8) 



1 -f (fc_Ro)"+3-^ 
where -Ro is a measure of the non-linear scale, defined from 
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0.55 



the linear variance aj^{Ro) = 1 with a Gaussian window 
function. The function Qzin) is given by 



Qain) 



(4-2") 
(l + 2"+i; 



(9) 



Note that at large scales, when the functions a = b = c = 1, 
we recover the tree-level FT expression for the bispectrum; 
on the other hand, at small scales a^ = (7/10) Q-i and 
6 = c = 0, so Eqs. (g|) and (|5|) give Q123 = Qsin) indepen- 
dent of triangle configuration: this is the prediction of hy- 
perextended perturbation theory (HEFT), which was shown 
to describe well the bispectrum (and higher order moments) 
in the strongly non-linear regime (SF). 

Summarizing, the non-linear bispectrum for scale-free 
intial conditions is given by Eq.(|2|) and Eqs.(Bl-H), together 
with the non-linear power spectrum obtained, for example, 
from the fitting formula of Peacock & Dodds (1996). In the 
next section we provide a generalization of the bispectrum 
fitting fornrula for CDM models. 



3 A FITTING FORMULA FOR CDM MODELS 

The numerical simulations used in this work correspond to 
cluster-normalized cold dark matter (CDM) models run by 
the Virgo collaboration (see e.g., Jenkins et al. 1998 for more 
details). There are four different models, SCDM, rCDM, 
OCDM, and ACDM, see table 1 for a description of the pa- 
rameters. These simulations each contain 256^ particles in a 
box 240^"^ Mpc on a side. They were run with an adaptive 
P'^M code (Couchman, Thomas & Pearce, 1995; Fearce & 
Couchman, 1997). There is only one realization for each cos- 
mological model, and all these have the same initial phases. 
We have also made measurements in a set of smaller simula- 
tions run by the Hydra Consortium (Couchman, Thomas & 
Pearce, 1995), with very similar results. The bispectrum was 
measured in this set of simulations as described in SCFFHM. 
Before generalizing the scale-free fitting formula (SF) 
reviewed in the previous section, a number of issues related 
to the CDM spectrum must be addressed. First, since the 
initial spectrum is no longer scale-free, there is no unique 
spectral index to be used in Eqs. (pI-h). Here we adopt the lo- 
cal spectral index obtained from the linear power spectrum. 
Another possibility would be to use at scale k the spec- 
tral index at scale fcz, obtained from the HKLM mapping, 
k = fei[l-|-47rfc^P(fc)]^/^ (Hamilton et al 1991), which maps a 
non- linear scale to a linear scale with the same number of en- 
closed pairs. Both approaches are reasonable but can only be 
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Figure 1. The bispectrum as a function of angle 6 between ki and k2 for scales as shown in the panels, for the OCDM model with 
0-8 = 0.85. 
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Figure 2. The bispectrum as a function of angle between ki and k2 for scales as shown in the panels, for the SCDM model with 
as = 0.61. 
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Figure 3. The bispectrum as a function of angle d between ki and k2 for scales as shown in the panels, for the ACDM model with 
0-8 = 0.55. 
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Figure 4. The bispectrum as a function of angle d between ki and k2 for scales as shown in the panels, for the rCDM model with 
0-8 = 0.61. 
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Figure 5. The bispectrum for equilateral triangles as a function of scale. Solid lines correspond to the fitting formula, symbols to the 
measurements in N-body simulations. 



simple approximations to a more complex situation; in prac- 
tice the transfer of power from large to small scales depends 
not only on the local spectral index but on the full shape 
of the power spectrum (e.g. Jain & Bertschinger 1994) and 
non-linear power at a given wavenumber is not just a func- 
tion of the linear power at another scale. In the transition to 
the non-linear regime, one-loop calculations (SCFFHM) and 
numerical simulations with initial power spectra cutoff be- 
yond some characteristic scale (Fry et al. 1993) suggest that 
the local index approach is a better choice. On the other 
hand, at small scales one would expect the HKLM choice 
to be clearly better, since for models such as CDM where 
the initial power spectrum decreases strongly with k, most 
of the power is generated by cascading from larger scales. 
Despite this, approximations for one-point higher order mo- 
ments in the non-linear regime based on the local spectral 
index such as EPT (Colombi et al. 1997) and HEPT (SF) 
show very good agreement even for CDM spectra. Further- 
more, attempts to use the HKLM choice gave poor results 
(Colombi et al. 1997). 

Another issue is the dependence of the fitting formula 
on cosmological parameters such as fim and Qa- Pertur- 
bative results show that the reduced bispectrum Q is very 
weakly dependent on cosmology at large scales, in particu- 
lar, F2(ki,k2) in Eq. rt3) is very well approximated by its 
Einstein-DeSitter value (Bouchet et al. 1995). In addition, 
SCFFHM showed that as long as ^^//^(nm) « 1, where 
f{Q,m) ~ f2m^, all higher-order moments become indepen- 
dent of cosmology at all orders in perturbation theory for 
a fixed erg- For these reasons, we shall assume that the fit- 
ting formula does not depend on cosmology; this is certainly 



a good approximation within the accuracy of our measure- 
ments in simulations. Thus, all the cosmology dependence 
of the bispectrum -6123 is encoded in the cosmology depen- 
dence of the non-linear power spectrum. 

We should also note that the asymptotic behavior of 
the bispectrum at small scales is not yet well tested against 
N-body simulations. Here we shall follow the standard as- 
sumption of the hierarchical ansatz, in which Q becomes 
constant in the deeply non-linear regime, as implicitly as- 
sumed in Eq. (0). This behavior is consistent with our mea- 
surements in simulations in the limited range of scales that 
we can test. 

To obtain a fitting formula for the bispectrum in CDM 
models we proceeded as follows. First, we used the fitting for- 
mula for the scale-free case by replacing kRo with g — k/kni , 
where 4TTk^iPL{kni) — 1, and Phik) is the linear power spec- 
trum at the corresponding redshift. This provides a reason- 
able, but not very accurate, fit to the simulation results, 
with an average deviation of davg = 40% from the N-body 
measurements. We then perturbed away from this fit by us- 
ing two methods. We introduced six free parameters, a^, 
such that q — > a^q in Eq. (m), n -|- 6 — > n -|- a2 in Eq. (H), 
0.2 -^ 0.2a3 in Eq. (Q), 4.5 -^ 4.5a4 in Eq. (|), q ^ a^q in 
Eq. d), and finally [0.7 Qsiri)]^''' -^ [0.7 Q3{n)Y/'^ap{z), 
which quantifies deviations from HEPT that depend on red- 
shift. We then obtained the best fit parameters by doing a x^ 
minimization between the predicted and measured bispec- 
trum for all triangles, using a large grid of models (contain- 
ing about 10^ of them). We also looked for the best fit using 
the Levenberg & Marquardt Method (Press et al. 1992), the 
solution was consistent with the grid approach but some- 



what sensitive to the initial guess for a^. As a result of this 
procedure we get a best fit model with parameters ai = 0.25, 
a2 = 3.5, as = 2, 04 = 1, as = 2, ae = —0.2, that is 



a{n, k) 
b{n, k) 
c{n, k) 



1 + cr: 



'{z)[0.7Qs{n)]'^'iq/4) 



n+3.5 



1 + (,?/4)"+3-5 

1 + 0.4 (n + 3) q"+'-^ 
1 + g"+3-5 

1 + 4.5/ [1.5 + (n + 3)''] (2g)"+3 



1 + (2g) 



n+3.5 



(10) 
(11) 
(12) 



and q = k/kni{z). The average deviation of this fitting for- 
mula from the measurements in the simulations is about 
rfavg = 15%. We also found similarly good fits with davg = 
16% for a\ = 0.25, a2 = 4, aa = 1.5, 04 = 1, as = 1.5 
and as = 0, so the small deviations from HEPT implied by 
Eq.(|lO[) do not seem statistically significant. 

Figures U-o illustrate the results. Although we have 
measured, and used in the fitting procedure, the bispectrum 
for many triangle shapes, we show for simplicity only tri- 
angles where ^2 = 2fci as a function of angle 9 between ki 
and k2, and equilateral triangles as a function of scale k. 
Since we only have a single realization for each cosmologi- 
cal model, we cannot show rigorous error bars, which in any 
case would only be part of the story, since difi'erent triangles 
in a single panel are correlated: highly so in the non-linear 
regime. However, we can obtain a reasonable estimate of the 
errors in a single realization by scaling the variance obtained 
from multiple realizations of ACDM lower resolution runs 
(Scoccimarro et al 2000b): this gives an error in Q of about 
±0.25 for scales k > 0.3 h/Mpc approximately independent 
of scale. 

Figures fll-W show a comparison of the fitting formula to 
the N-body measurements for triangles of difi^erent scales (as 
denoted in each panel) as a function of triangle shape. The 
reduced bispectrum Q shows the familiar configuration de- 
pendence at large scales (top-left panel), and then flattens 
out at small scales (bottom-right panel). This behavior is 
well reproduced by the fitting formula, and it is also well 
understood physically (SCFFHM): at large scales the bis- 
pectrum probes the shapes of large-scale structures, which 
leads to an anisotropic Q, whereas at small scales all config- 
urations become equally probable due to virialization. Note 
that the flattening is more evident in higher as models, as 
expected, since at a given scale these models have a stronger 
level of non-linearity. 

Comparison between corresponding panels in Fig. H and 
Fig. H and also between the top panels in Fig. H illustrate 
the lack of sensitivity of Q (for a fixed level of non-linearity, 
or as) to cosmological parameters such as Qrn and JIa, even 
in the non-linear regime. Note that in some cases the agree- 
ment at the smallest scales between the fitting formula and 
the simulations is not as good as at larger scales, especially 
in the rCDM case; however, this does not seem very serious 
given the estimated error of ±0.25 and that different trian- 
gles are strongly correlated. More numerical work would be 
necessary to assess this issue. 



4 CONCLUSIONS 

We provided a fitting formula for the non-linear evolution 
of the bispectrum in CDM models, which fits the measure- 
ments in N-body simulations to an accuracy of 15%. This 
will be a useful tool for making predictions about the non- 
linear behavior of non-Gaussian statistics related to the bis- 
pectrum of the density field, and is analagous to the role 
played by fitting formulas for the non-linear evolution of the 
power spectrum (Hamilton et al. 1991; Jain, Mo & White 
1995; Peacock & Dodds 1996). Indeed, Van Waerbeke et al. 
(2000) apply the results obtained here to the calculation of 
the skewness of the convergence field in weak gravitational 
lensing surveys. 

There are some issues in this work that merit further 
investigation. One obvious aspect is improving the accuracy 
of the fitting function; in view of the upcoming large galaxy 
and weak lensing surveys it would be useful to have a more 
accurate prediction for the non-linear evolution of the bis- 
pectrum. The main need in this regard is the availability 
of multiple independent numerical realizations of each cos- 
mological model, as the bispectrum is rather sensitive to 
the presence of massive clusters in the simulation volume 
(Colombi et al. 1996). Another issue is the extension to- 
wards smaller scales, and the validity of the hierarchical 
ansatz. Analytical models based on dark matter halos pre- 
dict that at small scales, k > 3 h/Mpc, Q should increase 
with k (Ma & Fry 2000, Scoccimarro et al. 2000b). Testing 
this prediction reliably will require the use of large- volume, 
higher-resolution simulations. 
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